function dxdt = dynamic_equation_modify(t, x, params)
    % params: 结构体，包含所有必要的参数
    % 从参数结构体中获取值
    m = params.m;          % 质量
    r_Gx = params.r_Gx;   % 惯性半径
    r_Gy = params.r_Gy;
    r_Gz = params.r_Gz;
    S_ref = params.S_ref; % 参考面积
    g_0 = params.g_0;     % 重力加速度
    rho = params.rho;     % 空气密度
    
    % 气动参数
    Cx0 = params.Cx0;
    Cy0 = params.Cy0;
    Cy_a = params.Cy_a;
    Cz0 = params.Cz0;
    Cz_b = params.Cz_b;
    
    % 力臂参数
    r_TBx = params.r_TBx;
    r_TBy = params.r_TBy;
    r_TBz = params.r_TBz;
    r_TDx = params.r_TDx;
    r_TDy = params.r_TDy;
    r_TDz = params.r_TDz;

    % 状态量分解
    rx = x(1); ry = x(2); rz = x(3);
    vx = x(4); vy = x(5); vz = x(6);
    q_0 = x(7); q_1 = x(8); q_2 = x(9); q_3 = x(10);
    omega_x = x(11); omega_y = x(12); omega_z = x(13);
    
    v = sqrt(vx^2 + vy^2 + vz^2);
    
    % 姿态矩阵计算
    C_BI = [1-2*(q_2^2+q_3^2), 2*(q_1*q_2+q_0*q_3), 2*(q_1*q_3-q_0*q_2);
            2*(q_1*q_2-q_0*q_3), 1-2*(q_1^2+q_3^2), 2*(q_2*q_3+q_0*q_1);
            2*(q_1*q_3+q_0*q_2), 2*(q_2*q_3-q_0*q_1), 1-2*(q_1^2+q_2^2)];
    
    % 速度转换到体轴系
    vb = C_BI * [vx;vy;vz];
    vbx = vb(1); vby = vb(2); vbz = vb(3);
    
    % 计算攻角和侧滑角
    alpha = atan2d(-vby, vbx);
    beta = asind(vbz / v);
    
    % 动压计算
    dp = 0.5 * rho * v^2;
    
    % 气动力计算
    Cx = Cx0;
    Cy = Cy0 + Cy_a * max(min(alpha,15),-15);
    Cz = Cz0 + Cz_b * max(min(beta,15),-15);
    
    Fxb = -Cx * dp * S_ref;
    Fyb = Cy * dp * S_ref;
    Fzb = -Cz * dp * S_ref;
    Fxd = 0; Fyd = 0; Fzd = 0;
    
    % 微分方程组
    dxdt = [vx;
            vy;
            vz;
            -((Fxb + Fxd)*(2*q_2^2 + 2*q_3^2 - 1) + (2*q_0*q_3 - 2*q_1*q_2)*(Fyb + Fyd) - (2*q_0*q_2 + 2*q_1*q_3)*(Fzb + Fzd))/m;
            - g_0 - ((Fyb + Fyd)*(2*q_1^2 + 2*q_3^2 - 1) - (2*q_0*q_3 + 2*q_1*q_2)*(Fxb + Fxd) + (2*q_0*q_1 - 2*q_2*q_3)*(Fzb + Fzd))/m;
            -((Fzb + Fzd)*(2*q_1^2 + 2*q_2^2 - 1) + (2*q_0*q_2 - 2*q_1*q_3)*(Fxb + Fxd) - (2*q_0*q_1 + 2*q_2*q_3)*(Fyb + Fyd))/m;
            - (omega_x*q_1)/2 - (omega_y*q_2)/2 - (omega_z*q_3)/2;
            (omega_x*q_0)/2 - (omega_y*q_3)/2 + (omega_z*q_2)/2;
            (omega_y*q_0)/2 + (omega_x*q_3)/2 - (omega_z*q_1)/2;
            (omega_y*q_1)/2 - (omega_x*q_2)/2 + (omega_z*q_0)/2;
            -(- m*omega_y*omega_z*r_Gy^2 + m*omega_y*omega_z*r_Gz^2 + Fyb*r_TBz + Fyd*r_TDz - Fzb*r_TBy - Fzd*r_TDy)/(m*r_Gx^2);
            (- m*omega_x*omega_z*r_Gx^2 + m*omega_x*omega_z*r_Gz^2 + Fxb*r_TBz + Fxd*r_TDz - Fzb*r_TBx - Fzd*r_TDx)/(m*r_Gy^2);
            -(- m*omega_x*omega_y*r_Gx^2 + m*omega_x*omega_y*r_Gy^2 + Fxb*r_TBy + Fxd*r_TDy - Fyb*r_TBx - Fyd*r_TDx)/(m*r_Gz^2)];
end